Source of Acquisition 
NASA Goddard Space Flight Center 


Post-Newtoiiia.ii Black-Hole Inspiral Initial Data witli Waves for Punctuie Simulations 

B. J. Kelly, 1,2 W. Tichy, 3 M. Campanelli, 4,2 and B. F. Whiting 5,2 

gravitational Astrophysics Laboratory , NASA Goddard Space Flight Center, 8800 Greenbelt Rd., Greenbelt, MD 20771, USA 
2 Center for Gravitational Wave Astronomy, Department of Physics and Astronomy, 

The University of Texas at Brownsville, Brownsville, Texas 78520 
3 Department of Physics, Florida Atlantic University, Boca Raton Florida 33431-0991 
i Center for Computational Relativity and Gravitation, 

School of Mathematical Sciences, Rochester Institute of Technology, 

78 Lomb Memorial Drive, Rochester, New York 14623 
5 Department of Physics, University of Florida, Gainsville Florida 32611-8440 

(Dated: March 9, 2007) 

We present improved post-Newtonian-inspired initial data for non-spinning black-hole binaries, 
suitable for numerical evolution with punctures. We revisit the work of Tichy et al. [1] and explicitly 
calculate remaining integral terms. Thereby we improve the accuracy in the far zone, by including 
realistic gravitational waves in the initial data. We investigate the behavior of this data both at the 
center of mass and in the far zone, demonstrating agreement of the transverse-traceless parts of the 
new metric with quadrupole-approximation waveforms. An advantage of these data is that they can 
be used for numerical evolutions to make a direct connection between the merger waveforms and 
the post-Newtonian inspiral waveforms. 
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I. INTRODUCTION 

Post-Newtonian (PN) methods have played a funda- 
mental role in our understanding of the astrophysical im- 
plications of Einstein’s theory of general relativity. Most 
importantly, they have been used to confirm that the ra- 
diation of gravitational waves accounts for energy loss in 
known binary pulsar configurations. They have also been 
used to create templates for the gravitational waves emit- 
ted from compact binaries which might be detected by 
ground-based gravitational wave observatories, such as 
LIGO [2, 3], and the NASA/ESA planned space-based 
mission, LISA [4, 5]. However, PN methods have not 
been extensively used to provide initial data for binary 
evolution in numerical relativity, nor, until recently (see 
[6]), have they been extensively studied so that their lim- 
itations could be well identified and the results of numer- 
ical relativity independently confirmed. 

Until the end of 2004, the field of numerical relativ- 
ity was struggling to compute even a single orbit for 
a black-hole binary (BHB). Although debate occurred 
on the advantages of one type of initial data over an- 
other, the primary focus within the numerical relativ- 
ity community was on code refinement which would lead 
to more stable evolution. Astrophysical realism was 
verjr much a secondary issue. However, this situation 
has radically changed in the last few years with the in- 
troduction of two essentially independent, but equally 
successful techniques: the generalized harmonic gauge 
(GHG) method developed by Pretorius [7] and the “mov- 
ing puncture” approach, independently developed by the 
UTB and NASA Goddard groups [8, 9]. Originally in- 
troduced by Brandt & Briigmann [10] in the context of 
initial data, the puncture method explicitly factored out 
the singular part of the metric. When used in numerical 


evolution in which the punctures remained fixed on the 
numerical grid, it resulted in distortions of the coordinate 
system and instabilities in the BSSN evolution scheme. 
The revolutionary idea behind the moving puncture ap- 
proach was precisely, not to factor out the singular part 
of the metric, but rather evolve it together with the reg- 
ular part, allowing the punctures to move freely across 
the grid with a suitable choice of the gauge. 

A golden age for numerical relativity is now emerg- 
ing, in which multiple groups are using different com- 
puter codes to evolve BHBs for several orbits before 
plunge and merger [11, 12, 13, 14, 15, 16, 17, 18]. Com- 
parison of the numerical results obtained from these 
various codes has begun [19, 20, 21], and compari- 
son with PN inspiral waveforms has also been carried 
out with tantalizing success [6, 22, 23]. The appli- 
cation of successful numerical relativity tools to study 
some important astrophysical properties (e.g. preces- 
sion, recoil, spin-orbit coupling, elliptical orbits, etc) 
of spinning and/or unequal mass-black hole systems is 
currently producing extremely interesting new results 
[24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38]. It 
now seems that the primary obstacle to further progress 
is simply one of computing power. In this new situation, 
it is perhaps time to return to the question of what initial 
data will best describe an astrophysical BHB. 

To date, the best-motivated description of pre-merger 
BHBs has been supplied by PN methods. We might ex- 
pect, then, that a PN-based approach would give us the 
most astrophysically correct initial data from which to 
run full numerical simulations. In practice, PN results 
are frequently obtained in a form ill-adapted to numeri- 
cal evolution. PN analysis deals with the full four-metric, 
in harmonic coordinates; numerical evolutions frequently 
use ADM-type coordinates, with a canonical decomposi- 
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tion of the four-metric into a spatial metric and extrinsic 
curvature. 

Fortunately, many PN results have been translated 
into the language of ADM by Ohta, Damour, Schafer 
and collaborators. Explicit results for 2.5PN BHB data 
in the near zone were given by Schafer [39] and Jara- 
nowski & Schafer (JS) [40], and these were implemented 
numerically by Tichy et al. [1] . Their insight was that the 
ADM-transverse-traceless (TT) gauge used by Schafer 
was well-adapted to the puncture approach. 

The initial data provided previously by Tichy et 
al. already includes PN information. It is accurate up 
to order (u/c) 5 in the near zone (r < A), but has the dis- 
advantage that the accuracy drops to order ( v/c ) 3 in the 
far zone ( r A) [here A is the gravitational wavelength] . 
Thus the data were incomplete in the sense that they did 
not include the correct TT radiative piece in the metric, 
and thus did not contain realistic gravitational waves. 

In this paper, we revisit the PN data problem in ADM- 
TT coordinates, with the aim of supplying Numerical 
Relativity with initial BHB data that extends as far 
as necessary, and contains realistic gravitational waves. 
To do this, we have evaluated the “missing pieces” of 
Schafer’s TT metric for the case of two non-spinning par- 
ticles. We have analyzed the near- and far-zone behavior 
of this data, and incorporated it numerically in the Cac- 
tus framework. 

The remainder of this paper is laid out as follows. In 
Section II, we summarize the results of Schafer (1985) 
[39], and Jaranowski & Schafer (1997) [40] and their ap- 
plication by Tichy et al. (2003) [l], to the production of 
puncture data for numerical evolution. In Section III, we 
derive briefly the additional terms necessary to complete 
h TT to order (u/c) 4 , deferring details to Appendix A. 
We pause to consider the issue of orbital phase calcula- 
tion in PN theory in Section IV. In Section V, we study 
the full data both analytically and numerically. Section 
VI summarizes our results, and lays the groundwork for 
numerical evolution of this data, to be presented in a 
subsequent article. 


incorporating an arbitrary number of spinless point par- 
ticles, with arbitrary masses For N particles, the 
lowest-order contribution to the conformal factor is: 

(i) 

z — 4 V A 

A~1 A 

where r A = y/x — x A is the distance from the field point 
to the location of particle A. In principle hj? is com- 
puted from 

h ij *ij ^ret S kh ( 2 ) 

where is the inverse d’Alembertian (with a “no- 
incoming-radiation” condition [42]), ski is a non-local 
source term and <5?™ is the TT-projection operator. 
In order to compute hj? we first rewrite it as 

K V = -^“[A- 1 + A"')] s u 

= A— (3) 

Note that the near-zone approximation h 4 j T ( KZ ) of h ] ] T 

has already been computed in [39] up to order 0(v/c) 4 
(see also Eq. 10 below). The last term in Eq. (3) is diffi- 
cult to compute because 

s ki = 167rG^ PAkPAI d(x - x A ) + 4^ (4) 

A m A 4 

is a non-local source. However, we can approximate sjy 

by 
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(5) 


II. ADM-TT GAUGE IN POST-NEWTONIAN 
DATA 

The “ADM-TT” gauge [39, 41] is a 3+1 split of data 
where the three-metric differs from conformal flatness 
precisely by a TT radiative part: 

9ij = + Vij+hf?, 

4 = 0. 

The fields d>, 7r y " and fiTT can all be expanded in a 
post-Newtonian series. Solving the constraint equations 
of 3+1 general relativity in this gauge, [39, 40] obtained 
explicit expressions up to 0(v/c) 5 valid in the near zone, 


and show that 

h ludiv) = ~ 5 7j Tkl ( D 7et~^~ 1 )(skl- Ski) ~ 0{v/cf (6) 

in the near zone. Furthermore, outside the near zone 

h iUdiv) ~ 1 / r ' 2 ’ so that hj^div) falls much faster than 
rest of hj?, which falls off like 1/r. Hence 

hi? = h]?™ - - A-‘)f H + hj? m , (7) 

where hj^ div ) can be neglected if we only keep terms up 
to 0(v/c) 4 . 

The full expression for fiTT for A T interacting point par- 
ticles from Eq. (4.3) of [39] is (making the replacement 
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1/16 7r — > G where necessary) 


The first term in (8), / l ^ T ( NZ ) 


can be expanded in v/c as 
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/i rr(NZ) = /i TT(4) +/i T T (5) +0(v / c )6. 


( 9 ) 


The leading order term at 0(u/c) 4 , is given explicitly by 
Eq. (A20) of [40]: 
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(10) 


where sab = r A + r B + r A B ■ The other two terms in 
Eq. (8) can be shown to be small in the near zone (r <C A, 
where the characteristic wavelength A ~ 100M for d ~ 
10M). However, hj^ (AZ) is only a valid approximation 
to h? 1 ' in the near zone, and becomes inaccurate when 
used further afield. 

Setting aside these far-field issues, Tichy et al. [1] ap- 
plied Schafer’s formulation, in the context of a black-hole 
binary system, and constructed initial data that are ac- 
curate up to 0{v/c ) 5 in the near zone. In particular, 
they noted that the ADM-TT decomposition was well- 
adapted to the use in numerical relativity of punctures 
[10] to handle black-hole singularities. 

The PN-based puncture data of Tichy et al. have not 
been used for numerical evolutions. This is in part be- 
cause these data, just like standard puncture data, do not 
contain realistic gravitational waves in the far zone. To 
illustrate this, we restrict to the case of two point sources, 
and compute the “plus” and “cross” polarizations of the 
near-zone approximation for hj?: 

h ? Z) = h^ NZ) eiel (ID 

h<?*> = h5 r(JVZ) 44- ( 12 ) 

For comparison, the corresponding polarizations of the 
quadrupole approximation for the gravitational-wave 


strain are given by (paraphrasing Eq. (3.4) of [43]): 

h+ = — — (1 + cos 2 9)(nMfG\v?^ 3 cos($gw), (13) 
r 

h x = cos 9(nMfGW? /3 sin($Giv)i (14) 
r 

where M — u 3 ^ 5 M is the “chirp mass” of the binary, 
given in terms of the total PN mass of the system M = 
mi + m 2 , and the symmetric mass ratio v — mi m 2 /M 2 . 
§Gw(t) and few (t) are the phase and frequency of the 
radiation at time t, exactly twice the orbital phase $(r) 
and orbital frequency fi(r)/27r, evaluated at the retarded 
time r = t r =t - r. The lowest-order PN prediction for 
radiation-reaction effects yields a simple inspiral of the 
binary over time, with orbital phasing given by 

■Hr) = «l,j - ie 5 /®, (15) 

n(T) = m 6 ' 3 ' 3 - <16) 

where 0 = v(t c - r)/5M, M and v are given below 
(14), and t c is a nominal “coalescence time”. Finally, the 
angle 9 above is the “inclination angle of orbital angular 
momentum to the line of sight toward the detector” ; that 
is, just the polar angle to the field point, when the binary 
moves in the x-y plane. 

Additionally, to evaluate (11-12), we need the trans- 
verse momentum p corresponding to the desired separa- 
tion r i 2 - The simplest expression for this is the classical 
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FIG. 1: Plus polarization of the quadrupole (black/solid) 
and near-zone (red/dashed) strains observed at field point 
r = 100M, 8 = 7t/4, 4> — 0. The binary orbits in the x-y 
plane, with initial separation is ri 2 = 10M, with a nominal 
coalescence time t c « 1180M. Both phase and amplitude of 
^tt (4) are ver y wron g outside the near zone. 


Keplerian relation, which we give parametrized by 


riz = M(MD)- 2 / 3 , 

(17) 

p = Mu(Mnf /3 . 

(18) 


In Fig. 1 we compare the plus polarization of the two 
waveforms (11) and (13) at a field point r = 100M, 
0 = 7r/4, <j> = 0, for a binary in the x-y plane, with 
initial separation ri 2 = 10M. The orbital frequency of 
the binary is related to the separation ri 2 and momenta 
p entering (11) by (17-18). In the lowest level of approx- 
imation, the binary has a nominal PN coalescence time 
t c « 1180M. As expected, both phase and amplitude of 

A. K are wrong outside the near zone. This means that 

TT ( 4 ) 

the data constructed from have the wrong wave 

content, but nevertheless these data are still accurate up 
to order («/-c) 3 in the far zone. 

It is evident from the present-time dependence of (10) 
that it cannot actually contain any of the past history of 
an inspiralling binary. We would expect that a correct 
“wave-like” contribution should depend rather on the re- 
tarded time of each contributing point source. It seems 
evident that the correct behavior must, in fact, be con- 
tained in the as-yet unevaluated parts of (8). The requi- 
site evaluation is what we undertake in the next section. 


III. COMPLETING THE EVALUATION OF hf? 

To move forward, we simplify (8) and (10) to the case 
of only two particles. Then, (8) reduces to: 
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( 20 ) 


where 


Hr A \u } 
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( 21 ) 


r 


Here, the “TT projection” is effected using the operator 
P- := 5\ — ki W jk 2 . For an arbitrary spatial vector u, 


lTT 


u c u d (P t c P“- 

u i U j o 


r Pa P cd ) 
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Details on the evaluation of these terms are presented 
in Appendix A. After calculation, we write the result as 
a sum of terms evaluated at the present field-point time 
t, the retarded time t r A defined by 


-2 


k 


t ~ t r A — r A (t r A ) =0, 


(23) 
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and integrals between t r A and t, 


where the three parts are given by: 


Hf t a P] 


-^tt A Pi A 
+-^TT A pi 


■ A [u\t r A ] 


*]» 


(24) 


Pi *] = — t — yr | [w 2 - 5 (u ■ n A ) 2 ] 5 l J + 2 u* tF + [3 (u ■ n A ) 2 - 5 -u 2 ] n) 4 

+12 (u • n^i) n A | , (25) 

H^ ta [u; t r A \ = — j[-2u 2 + 2 (u- n A ) 2 ] 5 lj + Au l u j + [2u 2 + 2(u- n A ) 2 ] n A n j A 
r A {t A ) L 

-8 (u ■ n A ) u ^ n j } \ , (26) 

H 1 t? t A [ u\t r A -* t) = -G ^ dr ^-^~ | [-5u 2 + 9 (u ■ fi A ) 2 ] 5 l j + 6u* u j - 12 (u ■ h A )u( z n A 
+ [9 u 2 — 15 (u ■ fij i) 2 ] n\ n\ j 

— G f dr ^ | [u 2 - 5 (u ■ h A ) 2 ] 5 zj +2 u 1 u j — 20 (u- h A )u ^ n A 

Jt r A r A {Tf l 

+ [-5 u 2 + 35 (u- n A ) 2 ] n A n j A j . (27) 


In Fig. 2, we show the retarded times calculated for 
each particle, as measured at points along the x axis, for 
the same orbit as in Fig. 1. We also show the corre- 
sponding retarded times for a binary in a stable circular 
orbit. Since the small-scale oscillatory effect of the finite 
orbital radius would be lost by the overall linear trend, 
we have multiplied by the orbital radius. 

A. Reconciling with Jaranowski 8z Schafer’s hf J T ^ 

From the derivation above it is clear that hj? includes 
retardation effects, so it will not depend solely on the 
present time. We might even expect that all “present- 
time” contributions should vanish individually, or should 
cancel out. In fact, it can be seen easily from (25) that 
the “f” part of the second and third terms of Eq. (20) 
exactly cancel out the “kinetic” part (first line) of Eq. 
(10). That is, we can simply remove that line to begin 
with, and use the “t r ” part instead. One may also inquire 
whether the “t” parts of the fourth and fifth terms of Eq. 
(20) above, 



FIG. 2: Retarded times for particles 1 and 2, as measured by 
observers along the x axis at the initial time t = 0, for the 
binary of Fig. 1. To highlight the oscillatory effect of the 
finite-radius orbit on t r , we first multiply by r, which is the 
average field distance. 
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similarly cancel the remaining, “potential” parts of Eq. 
(10). The answer is “not completely”; expanding in pow- 
ers of 1 /r, we find: 
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h TT (pot, 4) + h JT (pot, now) = & ™1™2 r l2 | ( g + M pj/2 _ 35 ^4) g f , _ 4 (j + 5 ^2) nm 

-5 (1 + 6 W 2 - 7 W 4 ) nu ny + 2 W (7 + 9 W 2 ) (n 12 < n v + n 12j n u )j + 0(l/r 4 )(29) 


where W = sin# cos {<j> - 3>(t)), and $>(t) is the orbital 
phase of particle 1. That is, the “new” contribution can- 

TT ('4') 

cels the 1/r and l/r 2 pieces of v ' entirely. More- 
over, the &TT term in Eq. (20) will be neglected from 
here on because it is small both in the near and the far 
zone [39]. 

We note here two general properties of the contribu- 
tions to the full hj? . 

1 . In the near zone h? T 4 ^ is the dominant term since 
all other terms arise from (Ore 4 — A ~ 1 )ski- Thus 
all other terms must cancel within the accuracy of 
the near-zone approximation. 

TT ('4') 

2. ' ’ is very wrong far from the sources; thus, 

TT (4) 

the new corrections should “cancel” K ’ entirely 
far from sources. Note, however, that while hij = 
-□“ e 4 Sfei depends only on retarded time, its TT- 
projection h,J? = 5j^ kl hki has a more complicated 
causal structure; E.g. the finite time integral comes 
from applying the TT- projection. [Proof: If we 
had a source given exactly by Ski, would 

depend only the present time, /iy would depend 
only on retarded time, and /iy T would (as we have 
computed) contain a finite time integral term.] 

Additionally, the full /iy T agrees well with quadrupole 
predictions, which we demonstrate in Section V. 

IV. HIGHER-POST-NEWTONIAN PHASE 
INFORMATION 

The calculations presented above lead to waveform am- 
plitudes that are accurate up to 0{v/c) A everywhere. 
However, the phase of the waveforms (via /i TT ) is critical 
for detection purposes, and to model this to high accu- 
racy, we desire that our initial-data wave content already 
encode the phase as accurately as possible. 

A. Circular or Inspiral? 

Waveform phase is a direct consequence of orbital 
phase. To lowest approximation, we can assume a binary 
moving in a closed circular orbit (for zero eccentricity). 
Although physically unrealistic - since radiation reaction 
will lead to inspiral and merger of the particles - this is 
perfectly consistent with our amplitude calculations. Up 
to 2PN order, we can have closed circular orbits, where 


the sources’ linear momenta p are related to the separa- 
tion ri 2 by, say, Eq. (2.4) of [1], However, a circular orbit 
has implications for the retarded radiation, since radia- 
tion observed at each point in the present time will have 
arisen from an orbit of the same size (with same momenta 
and accelerations). As more distant field points have ear- 
lier retarded times, the circular orbit approximation will 
get worse for those further field points. 

These relations are at the heart of the quadrupole ap- 
proximation for the observed wave strain (13-14). As the 
signal phase is so crucial for detector scientists, a popular 
compromise set of waveforms has been arrived at, called 
the “restricted post-Newtonian” approximation. In this 
approximation, we work with the quadrupole approxima- 
tion wave amplitudes (thus ignoring post-Newtonian am- 
plitude corrections), but feed these high-order-accurate 
PN phases. That is, we still use expressions (13) and (14), 
but with more complicated expressions for the phase 
$(r) + $(f c ) 1 

Expressions for the higher-PN orbital phase can be 
found in many references; for instance, Eq. (6.29) of 
[44] gives the orbital phase to high PN order, in radia- 
tive coordinates. Higher-order corrections to the phase 
in ADM-TT gauge can in principle be obtained through 
a contact transformation, or from direct ADM-TT calcu- 
lations [45]. 

Given this phase (and the implied orbital frequency Q), 
we can now derive the separation ri 2 and momentum p. 
For instance, in the manner of [46], we find to second PN 
(post leading) order: 


D2(fi) 

M 

= (Mft)~ 2 / 3 (3 V) 



(18 - 81i/ - 8p 2 ) rr\\2/3 

72 z (Afn) a/3 , 

(30) 

P(ty 

Mv 

= (MQ) 1/3 + ^ 15 6 " v \msi) 



+ (441 — 324,-^) (Mfi)! /, 

I A 

(31) 


B. Numerical Implementation of High-Order Phase 

For numerically constructing initial data, the primary 
input is the coordinate separation of the holes. This 


1 For consistency, one could similarly “upgrade” / gw , but this 
may not be necessary for signal phase determination. 



7 



FIG. 3: The xx component of the full for a binary with 
initial separation d = 10 M in a circular (black/solid) or in- 
spiralling (red/dashed) orbit. The orbital configuration is the 
same as for Fig. 1, apart from the Keplerian relations, where 
we have used the higher-order Keplerian relations (30-31). 
Note the frequency broadening at more distant field points. 


must be maintained exactly so that the punctures can 
be placed on the numerical grid. To ensure this, we in- 
vert Eq. (30) to obtain the exact Q r corresponding to 
our desired ri 2 - 

Now we use Eq. (16) with t = 0 to find the coalescence 
time t c that yields this Q r . Once we have obtained t c , 
we can find the orbital phase $ and frequency Q at any 
retarded time r directly from Eqs. (15-16), and the cor- 
responding separation r \2 and momentum p from Eqs. 
(30-31). 



FIG. 4: Plus and cross polarizations of the strain observed at 
field point r = 100M, 6 = 7t/4, (j> — 0. Both the quadrupole- 
approximation waveform (black) and the full (red) waveforms 
coming from hjj are shown. The orbital configuration is the 
same as for Fig. 1. 


results, as they should. 

After having confirmed that we have a PN three-metric 
Qij that is accurate up to errors of order 0(v/c) 5 , and 
that correctly approaches the quadrupole limit outside 
the near zone, we are now ready to construct initial data 
for numerical evolutions. In order to do so, we need the 
intrinsic curvature K^, which can be computed as in 
Tichy et al. [1] from the conjugate momentum. The dif- 
ference is that here we use the full hj? instead of the 
near zone approximation h[' r ^ to obtain the conjugate 
momentum [39]. The result is 


V. NUMERICAL RESULTS AND INVARIANTS 

In Fig. 3, we show a representative component of the 
retarded-time part of hf? for both circular and leading- 
order inspiral orbits. For both orbits, we use the ex- 
tended Keplerian relations (30) and (31); otherwise the 
orbital configuration is that of Fig.l. We can see that 
the cumulative wavelength error of the circular-orbit as- 
sumption becomes very large at large distances from the 
sources. This demonstrates that using inspiral orbits 
instead of circular orbits will significantly enhance the 
phase accuracy of the initial data, even though circular 
orbits are in principle sufficient when we include terms 
only up to 0(v/c) 4 as done in this work. From now on 
we use onlj'' inspiral orbits. 

Next, we compare our full waveform hj^ (expressed 
as the combinations h + and h x ) at an intermediate-field 
position (r = 100M, 9 = 7r/4, 0 = 0) to the lowest-order 
quadrupole result. In Fig. 4, the orbital configuration is 
the same as for Fig. 1. As one can see, both the + and x 
polarizations of our hj? agree very well with quadrupole 


K ij 


-'•I’pn 


7T, 


(3) 


+ + (^(2)^(3) ) TT 


+0(v/c) 6 , 


(32) 


where the error term comes from neglecting terms like 
hj^ div) at 0{y/cf in and where -ip PN , and 
0(2) can be found in Tichy et al. [1]. An additional dif- 
ference is that the time derivative of is evaluated 
numerically in this work. Note that the results for gjj 
are accurate up to 0(v/c) 4 , while the results for Kij are 
accurate up 0(v/c ) 5 , because Kjj contains an additional 
time derivative [1, 47, 48]. 

Next we show the violations of the Hamiltonian and 
momentum constraints computed from gij and K,j , as 
functions of the binary separation ri 2 . As we can see in 
both panels of Fig. 5, the constraints become smaller for 
larger separations, because the post-Newtonian approx- 
imation gets better. Note that, as in [1], the constraint 
violation remains finite everywhere, and is largest near 
each black hole. 





FIG. 5: Upper panel: Hamiltonian constraint violation along 
the y axis of our new data in the near zone, as a function of 
binary separation ri 2 - Lower panel: Momentum constraint 
(y-component) violation of the same data along the x axis. 
The orbital configuration is that of Fig. 3. 


A. Curvature Invariants and Asymptotic Flatness 

In analysis of both initial and evolved data, it is often 
instructive to investigate the behavior of scalar curva- 
ture invariants, as these give some idea of the far-field 
properties of our solution. We expect, for an asymptoti- 
cally flat space-time, that in the far field, the speciality 
index S = TU 2 /! 3 will be close to unity. This can 
be seen from the following arguments. Let us choose a 
tetrad such that the Weyl tensor components ipi and ip 3 
are both zero. Further, we assume that in the far field 
ip 0 and ipi are both perturbations of order e off a Kerr 
background. Then 

<S«l-3^|i + 0(e 3 ), (33) 

Y2 

which is indeed close to one. Note however, that this 
argument only works if the components of the. Weyl ten- 
sor obey the peeling theorem, such that 1 P 2 ~ 0(r~ 3 ), 
ipo rsy 0(r 5 ) and ip a ~ 0(r J ). In particular, if ipo falls 
off more slowly than 0(r -5 ), S will grow for large r. 
Now observe that ipo ~ 0(r“ 5 ) rv-' M 3 /r 5 is formally of 
0(u/c) 6 . Thus in order to see the expected behavior of 
S « 1 in the far-field we need to go to 0(v/c)°. If we 
only go to 0(v/c) 4 (as done in this work) ipo consists of 
uncontrolled remainders only, which should in principle 
be dropped. When we numerically compute S we find 
that for our data, S deviates further and further from 
unity for large distances from the binary. This reflects 
the fact that the so-called “incoming” Weyl scalar ipo 


only falls off as 1 /r 3 , due to uncontrolled remainders at 
0 (u/c) 6 , which arise from a mixing of background and 
TT waveform. 


VI. DISCUSSION AND FUTURE WORK 

The exploration and validation of PN inspiral wave- 
forms is of crucial importance for gravitational-wave de- 
tection and our theoretical understanding of black-hole- 
binary systems. The goal of this work is to provide a 
step forward in such understanding by building a direct 
interface between PN approximations and numerical evo- 
lutions along the lines already developed in Ref. [1]. In 
this paper we have essentially completed the calculation 
of the transverse-traceless part of the ADM-TT metric 
to 0(v/c) 4 provided in [ 1 ], yielding data that on the ini- 
tial Cauchy slice will describe the space-time into the 
far-field. We have incorporated this solution into a nu- 
merical initial data routine, adapted to the “puncture” 
topology that has been so successful recently. We have 
demonstrated this data’s numerical properties on the ini- 
tial slice. 

The next step is to evolve this data with moving punc- 
tures, and investigate how the explicit incorporation of 
post-Newtonian waveforms in the initial data affects both 
the ensuing slow binary inspiral of the sources and the 
“new” radiation released from the system. We note es- 
pecially that our data are non-eonformally flat beyond 
0(u/c) 3 . We expect our data to incorporate smaller 
unphysical initial distortions in the black holes than is 
possible with conformal flatness, and hence less spurious 
gravitational radiation during the numerical evolution. 
We see this as a very positive step toward providing fur- 
ther validation of numerical relativity results for multiple 
orbit simulations by comparison with PN results where 
these are expected to be reliable. These initial data will 
also allow to fully evaluate the validity of PN results for 
merging binaries comparing them with the most accurate 
numerical relativity results. 
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APPENDIX A: DETAILS OF INTEGRAL 
CALCULATION 

Here we present some more details of the calculations 
that lead to the three contributions to Eq. (21): Eqs. 
(25-27). Inserting Eq. (22) in the general integral (21), 
we can write H^ TA [u] as a combination of scalar and 
tensor terms: 


H? TA [u] = 16t tG 


J dr | Ui Uj 


■ Sa 


Ia 


+ 


hj a + 


u c Ud 


red t . 
1 A 3 


Uc U d 


- [2 M c U(,;] t Jj) A + 
where the “I” integrals are defined as: 

r _ f d 3 kdw ( uj/k ) 2 e ikrA c ° s6 ~ i “ T 

lA = J 72 ^ k 2 -(w + ie) 2 

ij __ r d 3 kdu> k l W 

A = J (2tt) 4 k 2 




(to/k) 2 e 


2 0 ikrA cos$—iojT 


rijed 
1 A ~ 


/ 


k 2 — (u> + i e) 2 
d 3 /c dio k l y k c k d 


(27r) 4 fc 4 
(u/k)^ e^^ rA cos ®~i° J T 
k 2 — (uj + i e) 2 


(A2) 


(A3) 


(A4) 


Here T = t — r, and = x — xa- We have also taken 
our integration coordinates such that ta lies in the 2 di- 
rection, so that the dummy momentum vector k satisfies 

k-rA=krACOs8 , d 3 k = k 2 dk sin 9 d9 dp. 

Define the unit orthogonal vectors ft a = (0, 0, l) , £ = 
(cos <f > , sin (j), 0). Then we can write 

k = k cos 9 Ua + k sin 9 1 => k ■ ?a = Dt & ■ Ua- 
We can also define a projector tensor onto i\ 

Qab = S ab -n a n b => Q a b = 6 a b -n a n b 

=> Q a c Q c b = Q a b , Q a b n b = o , Q a b i b = r. 


1. Angular integration 

We will neglect the A subscript for now, until it be- 
comes relevant again. To calculate the integrals (A2-A4), 
we begin with the p integration. The only <j> dependence 
comes from the l parts of the k terms. It can be seen 


from elementary trigonometric integrals that: 

o, 

J dp t £ b = ixQ ah , 

J dpt l b l c l d = J (Q ab Q cd +Q ac Q bd + Q ad Q bc ) 
We 

I A bcd . Define w = cos 9. Then 


We use these to calculate the <j> integrals for I A b and 


d(f> 


dtj) 1 = 2 7T, 

k a k b 


k 2 


2 7r w 2 n a n b + 7r (l — w 2 )Q a b , 


' b a b b h c h d 

dtp — j = 2 7T w A n a n b n c n d 

K 


+6 7r w 2 (1 — w 2 ) Q^ a b n c n ^ 

+ ^(1 - w 2 ) 2 Q( ah Q cd \ 

So the next integrals will differ in their 9 dependence, 
contained in the powers of w above. The 9 integrals will 
contain the following basic types: 


»(«) s / +1 2 2HL 5 

7-i a 

r + 1 


(A5) 


92(a) 


94(a) 


dw w e' 


2 a w 


sinh a , cosh a 


-4- 


/-1 


a 


+4: 


sinh a 


(A6) 


r+1 , 4 aw „ sinh a cosha 

dww e =2 t 


a 


a 


2 ’ 


„ sinha cosha , sinha , . 

+24 — 5 48 — 2 b 48 — (A7) 

a 6 a 4 a b 

Now I ab and I abcd can be written as the linear combi- 
nations: 


a b 


r 


rabed 


1 


[Q ab l] r + 


(n“ n b — -Q ab )K 


(A8) 


i a n b n c n d -3Q (ab n c n d) + ^Q (ab Q cd A L 


+ 


+- 


3 Q {ab n c n d) - | Q (ab Q cd) ) K 


Q(abQCd) J 


(A9) 


I here can be expressed in terms of go (a) above: 


(ko 


d 3 k 


;/kf 


_?■ k r cos 0 — iuj T 


(2 7 r) 3 J 2 7r k 2 — (a; -f i e) 2 

^ 2 ^ r d k 


(2 7r)’ 
du) 
(2tt) : 


■ c 0 e 


1/2 


fc 2 — (u + i e) 2 


g 0 (ikr) 


J2 e~ iuT j Q _ (A10) 
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The 1 /2 factor is because we moved to integrating k over 
the whole real line instead of the positive half-line (this 
is permissible as g n (a) is an even function of a). K and L 
are defined analogously to 7, but with extra even powers 
of cos 6 = w: 


K = 


L = 


du> 


/ 


d 3 k (bj/k) 2 e ikr cose-iuiT CQS 2 t 


(2tt) 3 J 2 7T k 2 — (u) + i e) 2 


,2 -iwT r °° 
u> e 


1/2 


(2tt) : 

dbJ 


I - oo dk k 2 -(u> + i e) 2 


u> 2 e ~ iu,T 


•72, 


(2 7r) 3 

du fd 3 k (u/k) 2 
(2tt) s 


r d'‘k 0 

J 2 7T k 2 — 


(w -F i e) 2 


g 2 .{ikr) 
(All) 

ikr cos 0—iioT ^^4 ^ 


^ w 2 e~ iwT r °° 


1/2 


(2 7r) 3 
dhJ 

(2 7r) 3 


I- oo * k 2 - (u + i e) 


w 2 e~ iwT 


J 4 . 


2 54(ifer) 

(A12) 



2. Momentum integration 


Now we address the k integrals, defined as: 


FIG. 6: Contours needed to complete integration over k (top) 
and u> (bottom). 


dk 


fn(k) = f 


dk fn (k) 


dkf n (k), 


the J n : 


where we collect the positive exponents in the g n in the 
integrand of /+ (k), and the negative exponents in f~ ( k ): 


f+m = 9%(ikr)/2 

JnW- + 


fn(k) = 


g n (ikr)/2 
k 2 — (uj + ie) 2 ' 


We calculate this as the sum of contour integrals of the 
“plus” and “minus” integrands (necessary, as the oppo- 
site signs require different contours). Each of these has 
poles at k = 0, k = k + = u + i t, and k = fc_ = — u> — it 
(the first of these is from the g n ). We integrate the “plus” 
integrands anticlockwise around the contour C\ (blue), 
and the “minus” integrands anticlockwise around the 
contour C 2 (red) (see Fig. 6); taking the limit |fc| — > 00 , 
the contribution from the curved segments vanishes, and 
the residue theorem gives us: 


7 0 

7 2 

Ja 


We ir(w+i€) 7 

r(u) + ie) 2 r(u> + i e) 


2 5 


(A14) 


(A15) 


ne ir (w+ie) Keirjio+ie) [-2 + 2 i r (w + i e)] 

r(u> + i e) 2 r 3 (u> + i e) 4 

2 7T 

+r 3 (uj + i e) 4 ’ 

ne ir(v+ie) 4 7re ir ( w + ie ) 

r(u) + ie) 2 r 5 (to + ie ) 6 

-3 r 2 (w + i e) 2 + ir 3 (u> + i e) 3 ] 

24 7T 


r 5 (w + i e) 6 ’ 


[6 — 6 i r (u> + i e) 

(A16) 


Jn = 2 xi Res [/+, k + ] - 2 7 ri Res[/„ , k-] 

+7riRes[/,J',0] — 7riRes[/“,0]. (A13) 


Calculating the residues, we find the values of each of 


3. Frequency integration 


Now we perform the w integration. Inserting the re- 
sults (A14-A16) into (A10-A12) respectively, we see that 
each of 7, K and L contains a delta function, which we 
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can extract: 


I = 
K = 


L = 


1 


4nr 

1 

4nr 

r 

+ 


[S(T-r)-6(T)), 


5(T — r) + e 


doj 

(2 7r) 3 


dw 

J (2 5r) 3 

iu)T F 2b (u), 

dw 


-iu> (T — r) 


F 2a (uj) 


4?rr 


<5(T-r)+e— J e -*- (T-r) 


+ 


dw 

(2^ 


—ioj T 


Fa fc(w), 


where the new terms on the RHS come from the Js above, 
grouped by exponential, as that is what determines the 
contours chosen during integration (see Fig. 6): 


The only pole is at w = —ie, so if we can close in the 
upper half-plane, we’ll get zero. 


®.For T < 0, both the “a” and “b” integrals can be 
closed in C\. Result; zero contribution. 


® For 0 < T < r, the “a” integrals can be closed 
in Ci, but the “b” integrals must be closed in C 2 . 
Result: “b” contribution. 


• For T > r, both the “a” and “b” integrals must be 
closed in C 2 . But then the “a” and “b” residues 
cancel out. Result: zero contribution. 


F2a{u) 

F2b(u) 

^4a(w) 


Fib(uj) 


7r w 2 [—2 + 2 i r (w + i e)[ 
r 3 (w + i e) 4 
2nu> 2 


r 3 (w + ie) 4 ’ 


7 TUI 


[24 — 24 ir(u> + i e) 


r 5 (w + i e) 6 
— 12 r 2 (w + i e) 2 + 4 i r 3 (w + i e) 3 ] , 
247TW 2 
r 5 (u + % e) 6 


Now the residues are as follows (taking the e — » 0 limit): 


Res 


Res 


[e-^(T- r )F 2a (u;),_ieJ 
Res [e"* w T F 2 b(to), - i e] 
e ~ iuJ (T- r ) F i a (u), —i e 
Res [e~ iwT F 4 (,(w), — ie] 


2n iT 
pi 1 
2iriT 

’ 

47riT 3 

^5 ’ 

47riT 3 


Thus the only interesting contribution happens in the 
interval 0 < T < r t — r(r) < r < t. In this case, the 
final integrals yield 


dw 

(2 7r) 3 
du> 


-iu> ( T—r ) 


ioj (T—r) 


F 2b {u) 


Fibiu) 


J ( 2 tt ) 3 

leading to the final result for K and L: 


T 

2 7rr 3 ’ 

T 3 


I = 
K = 
L = 


- S(T-r ) 


S(T), 


4nr~' y ’ ' 4nr 

—8 (T-r)-Q(T)Q(r-T) T 
4nr 

1 


47 rr 


2 nr 3 ' 

8(T — r) — Q(T) Q(r — T) 


We use these to calculate the Pi and pi kl : 


pi = 

jij kl 


rd n° 


— ^ — 5(T — r) — @(T) 0(r — T) T 
An r 




(A17) 


n ! n k n l 


— — <5(T — r) — 0(T) Q(r — T) 
4nr nr 


R) -3Q«>n k n‘^(T)@(r-T) (W 


r 3 


n r° 


1 j(T) + e(T)0(r-T) 

8 \ 4 7rr 


T 


7rr° 


jp3 

nr 5 


(A18) 


4. Time integration 


The final integrations will be over the source time r. The “crossing times” for the two 0 functions are r = t and 
r = t r , where t is the present field time, and t r the corresponding retarded time defined by (23). Now taking a general 
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function y(r), we find that 


/ OO 

dr I A v(t) = 

-OO 


-OO 

fOO 


dr I a y(r) 


dr I 1 2 / kl y(r) 


y(t r A ) y(t) 

4irr A (t r A ) 4Trr A (t)’ 




„k „l y( T ) 

n A n A n A n A 4 n Ta 


A 4 nr A 


J T=t T . 


T=t Jt A 

n kl ) y ( r ) 
4 irr A 


J dr ^3 n l A n\ — 


ij\ (' t-r)y(T ) 

" 1 47rr A (r) 3 4 ’ 


QTQa 


r=t 


+ 


+ 


l^dT ( 3 Q { A n An l A + 4 Q ( aQa 1) ') 2 ir A {l) 3y ^ 

J tr dr (^-n\n\n k A n l A + 3Q { A : n k A n l l-^Q ( A ° Q A l) ^ ^jy(r). 


These can now be substituted into the general integral (Al). We write the result as a sum of terms at the present 
field-point time t, the retarded time t r A , and interval terms between them, 

4t ^ [u] = A [tf; t] + A [n; t r A ) + H & A [S; t r A -*t], 


Htt a \ u 'i A ~~ 


4 G 

r A (t) 

2 u k u c * 


u* v? - — 8' 3 
1 0 j)fe 


+ 


5« + 


UkUl 








4G 

r A{t r A ) 


u l v? <W 

2 


2 Uk n 3 ) n A + 

J t\ 


+ 

1 1 

y] + [ 

J t A 

-u k Uf 

L 2 . 

r n l A n A n A n A ] 
t7 A J 


Ufa Ul 


n A n\ 


-^tt a\ u ' ^ 


t] = -4 G f* dr ^ 
Jt r . 


r A (r) 3 


2 Mfe 


3 n A n A — K j + 


(3 n A n k A 
(t~r) s 


3ri A n A — S lj ) + 
Uk ui 


' U k Ui 
- 2 


\n A n 


a 11 A 


(3- fc - l 

(ij „k l) ^ 


S kl ) 8 ij 


6 Q a n A n A - ~ Q a Q a 


' UG L dT rA(r)S 


^J 1 ] ( n\ n\n\- 3 Q ( A 3 n k A n l) A + - Q ( A 3 Q A l) 
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